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Abstract 

Maerl beds are sensitive biogenic liabitats built by an accunnulation of loose-lying, non-genlculate coralline algae. While 
these habitats are considered hot-spots of marine biodiversity, the number and distribution of maerl-forming species is 
uncertain because homoplasy and plasticity of morphological characters are common. As a result, species discrimination 
based on morphological features is notoriously challenging, making these coralline algae the ideal candidates for a DNA 
barcoding study. Here, mitochondrial (C0I-5P DNA barcode fragment) and plastidial (psfaA gene) sequence data were used 
in a two-step approach to delimit species in 224 collections of maerl sampled from Svalbard (78°96'N) to the Canary Islands 
(28°64'N) that represented 10 morphospecies from four genera and two families. First, the C0I-5P dataset was analyzed with 
two methods based on distinct criteria (ABGD and GMYC) to delineate 16 primary species hypotheses (PSHs) arranged into 
four major lineages. Second, chloroplast (psbk) sequence data served to consolidate these PSHs into 13 secondary species 
hypotheses (SSHs) that showed biologically plausible ranges. Using several lines of evidence (e.g. morphological characters, 
known species distributions, sequences from type and topotype material), six SSHs were assigned to available species 
names that included the geographically widespread Phymatolithon calcareum, Lithothamnion corallioides, and L glaciale; 
possible identities of other SSHs are discussed. Concordance between SSHs and morphospecies was minimal, highlighting 
the convenience of DNA barcoding for an accurate identification of maerl specimens. Our survey indicated that a majority of 
maerl forming species have small distribution ranges and revealed a gradual replacement of species with latitude. 
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Introduction 

Maerl or rhodolith beds are accumulations of slow-growing, 
unattached non-geniculate (non-articulated) coralline algae that 
buUd three-dimensional habitats [1] that accommodate a wide 
biodiversity and are, therefore, considered as hotspots of marine 
life [2]. Commercial dredging together with a range of indirect 
impacts (bottom-fishing, aquaculture, eutrophication, sediment 
dredging) are known to negatively afTect their conservation and 
structure [3] . As a result, maerl beds are listed as threatened and/ 
or declining habitats by OSPAR (The Convention for the 
Protection of the marine Environment of the North-East Atlantic) 
[4] and treated as Special Areas of Conservation by EU Habitats 
Directive (Annex I, categories "sandbank covered by seawater all 



the time" and "large shallow inlets and bays"). In addition, the two 
coralline algal species commonly regarded as the main constituents 
of maerl beds in Europe {Phymatolithon calcareum (Pallas) W.H. 
Adey & D.L. McKibbin and Lithothamnion corallioides (P.L. & 
H.M. Crouan) P.L. & H.M. Crouan) are listed in Annex V as 
species whose eventual exploitation must be compatible with 
maintaining a favorable conservation status. 

Maerl beds are widely distributed along the coasts of the North- 
East Atiantic protected by the OSPAR Convention (OSPAR 
maritime area) and the adjacent Macaronesia. They are partic- 
ularly frequent in Scotiand, Ireland, Brittany and Galicia [4] at 
depths ranging from the intertidal to 50 m, but they reach up to 
60 min the Canary Islands and Madeira [4—7]. Up to 24 species of 
maerl have been recorded along the OSPAR area and southern 
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adjacent regions (Madeira Archipelago and Canary Islands) 
(Table 1). However, the actual number might be smaller as the 
taxonomic validity of eight taxa seems dubious because they were 
only reported in pioneer works from the 19* century and early 
20* {Lithothamnion hreviaxe Foslie, L. fornicalum Foslie, L. 
fruliriilosum (Ktitzing) Foslie, L. intermedium Kjellman, L. 
nodulosum Foslie, L. norvegicum (Areschoug) Kjellman, L. 
tusterense Foslie and L. ungeri Kjellman). Maerl-forming algae 
belong to six genera (Lithothamnion, Lithophyllum, Mesophyllum, 
Neogoniolilhon, Phymalholilhon and Spongiles) from two families 
(Corallinaceae and Hapalidiaceae) within the order Corallinales 
(Rhodophyta). According to the literature, four species are widely 
distributed and seemingly foUow a latitudinal replacement cline: 
Lithothamnion tophiforme (Esper) Unger and Lithothamnion 
glaciale Kjellman are mostly arctic and subarctic species, while 
P. calcareum and L. corallioides occur from the North and Celtic 
Seas to Madeira-Canary Islands {L. corallioides) or Azores (P. 
calcareum). The remaining 12 species occupy narrower latitudinal 
ranges. Eight of them were reported for regions with a long 
tradition of taxonomic surveys: Scodand, Britain, Ireland, and 
French Brittany [Lilholhamnion lemoineae Adey, L. sonderi 
Hauck, Phymatoliihon purpureum (P.L. Crouan & H.M. Crouan) 
Woelkerling & L.M. Ir\inc\ Mesophyllum lichenoides (J. EUis) Me. 
Lemoine, Lithophyllum dentatum (Kutzing) Foslie, L. duckerae 
Woelkerling, L. fasciculatum (Lamarck) Foslie and L. hibernicum 
Foslie). The remaining four include species reported for Macar- 
onesia [Neogoniolithon brassica-florida (Harvey) SetcheU & L.R. 
Mason, Lithophyllum crouanii Foshe, Spongiles fruliculosa Ktitz- 
ing) plus Mesophyllum sphaericum V. Peila, Barbara, W.H. Adey, 
Riosmena-Rodriguez & H.G. Choi, a maerl alga known from a 
single location in Galicia. An overwhelming majority of the 
previous studies have entirely relied on traditional practices of 
taxonomy based on morphological/ anatomical characters even 
though morphological identification of non-articulated coralline 
algae is challenging because phenotypic plasticity and convergence 
have resulted in a lack of well-defined diagnostic characters [8] . 
Only very recendy, DNA information has been used to identify 
and delineate European maerl-forming species [9] shedding light 
on our fragmentary knowledge on alpha diversity and genuine 
distribution of maerl-forming .species. 

The onset of the 2r' century has witnessed notable technolog- 
ical advances that can facilitate and accelerate the description of 
biodiversity [10,11]. In particular, DNA barcoding (http://www. 
ibol.org/) employs short, standardized DNA fragments as a 
diagnostic tool for identifying species [12]. In Rhodophyta, 
DNA barcodes obtained by sequencing the ,5' end of the 
mitochondrial gene cytochrome oxidase I (COI-5P) [13,14] 
proved very effective to shortcut the difficulties of morphology- 
based identification, allowing an accurate identification of known 
species [L5-23] and/ or the detection of cryptic ones [16,20,24— 
26]. In c:omparison, COI-5P sequences have been less frequentiy 
used to delineate new species of red algae [15,18,20,22,25]. 
Indeed, when DNA barcoding suggested the existence of new 
species, it was rarely regarded as a definitive proof; instead, it was 
used along with other genetic, morphological, geographical or 
ecological features in what has been referred to as integrative 
taxonomy [10,27,28]. 

Despite the above, DNA barcodes have been used as an 
exploratory tool for poorly surveyed taxa provided that the groups 
delineated by barcodes are regarded as primary species hypothesis 
(PSHs) [29,30]. PSHs can then be further tested with other sources 
of molecular, morphological, geographical and/or ecological 
evidence and even a multistep approach has been proposed to 
turn PSHs into more conclusive secondary species hypotheses 



(SSHs) [10] (for a similar approach see the molecular-assisted 
alpha taxonomy in [26]). In this context, the initial step is crucial 
and consists of the partition of COI-5P sequences into a set of 
PSHs. Recently, two methods based on distinct criteria have been 
proposed to infer the limits of the various PSHs when only 
molecular data are available and with no need for prior 
assumptions. On the one hand, the Automatic Barcode Gap 
Discovery (ABGD) [30] is a fast method that uses distances to split 
the sequence alignment into a set of PSHs following a recursive 
procedure until there is no further partitioning. This procc'dure 
automatically finds breaks in the distribution of genetic pairwise 
distances, referred to as the 'barcode gap', even when intra- and 
interspecific distances overlap. On the other hand, the General 
Mixed Yule Coalescent (GMYC) model [29] is based on detecting 
the shift of the branching rate that takes place in clock-constrained 
calibrated trees at the point of transition from species-level 
(speciation) to population-level (coalescence) evolutionary process- 
es. Using a likelihood criterion, the GMYC method permits an 
automated species delineation with appropriate statistical mea- 
sures of confidence. A later extension of the method allows for a 
variable transition from coalescent to speciation among lineages 
[31]. GMYC has been shown to be robust to a range of departures 
from its assumptions (varying population sizes among species, 
alternative scenarios for speciation/extinction, population growth 
and subdivision within species) but the accuracy of its delimitations 
can be compromised in groups with large effective population sizes 
and short divergence times between species [32] . Other potential 
shortcomings of the GMYC method have been extensively 
discussed elsewhere [33]. 

In this study, COI-5P sequences were obtained for maerl- 
forming species along the OSPAR maritime area and the adjacent 
Macaronesia. DNA barcodes were used to delimit a set of PSHs 
that were subsequently corroborated or challenged with indepen- 
dent molecular, geographic, and morpho-anatomical evidence. 

Materials and Methods 

Study area and sample collection 

As the study did not involve endangered or protected species, no 
specific permissions were required for sampling at most locations 
(see Table SI for coordinates). Still, sampling at two locations 
situat(;d within a national park in NW Spain (lat 42.211° long 
-8.896° and lat 42.394° long -8.815°) was conducted witii tiie 
permission of the park authority (Parrjue Nacional Maritimo 
Terrestre de las Islas Atlanticas de Galitda) and the park authority 
has signed a document stating its interest in the results of this 
study. 

Collection information for all the specimens used in this study is 
available at the Barcode of life Data Systems (BOLD: www. 

boldsystems.org; project "maerl-NE Atlantic"). From 1999 to 
2011, maerl specimens were extensively sampled by SCUBA 
diving or dredging within 4 out of the 5 regions of the OSPAR 
maritime area (Table SI); sampling ranged from the low intertidal 
to 40 m depth. Despite our efforts, no sample could be obtained 
for region V where maerl beds are probably restricted to the 
Azores Archipelago. To circumvent this shortage, samples were 
collected from the other two Macaronesian Archipelagos: the 
Canaries and Madeira. Sampling sites included type/neotype 
localities for 3 out of the 4 widely distributed maerl-forming 
species: L. corallioides (Rade de Brest, Finistere, France [34]), P. 
calcareum (Falmouth Harbour, Cornwall, England [35]), and L. 
glaciale (Spitsbergen Island, Svalbard Archipelago [36]). Addi- 
tionally, our samples included holotype material of the recentiy 
described M. sphaericum from the herbarium SANT of Uni- 
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versidad de Santiago de Compostela [9] and neotype material of 
P. calcareum from the herbarium BM of the British Museum of 
Natural History [35]. 

Freshly collected material was transported to the laboratory in 
seawater, oven-dried or air-dried as soon as possible, and 
vouchered in silica. Vouchers were temporarily deposited in a 
personal collection (BioCost Research Group, University of A 
Corufia, Spain) and will be transferred to the official SANT 
Herbarium. When feasible, several specimens per morphotype 
(differences in size, shape, branch thickness, and general habit) 
were sequenced at each locality. This sampling regime was 
intended to maximize the detection of species encountered at the 
various collecting sites while keeping the sequencing effort to a 
reasonable size. All specimens were photographed and identified 
to the lowest taxonomic level possible using morphology-based 
keys and specialized literature. A selection of specimens was also 
examined under scanning electron microscope (SEM, model 
JEOLJSM 6400, University of A Coruna). 

Field identification 

A putative species name was assigned to all specimens based on 
their gross morphology. Our 224 collections were partitioned into 
10 different morphospecies belonging to 4 genera and 2 families 
(Hapalidiaceae and Corallinaceae). Most plants were identified 
either as P. calcareum (140 collections) or as L. corallioides (60), 
two main constituents of maerl in Atlantic Europe. A much 
smaller number of collections fitted the description of L. glaciale 
(9), M. sphaericum (4), L. dentatum (2), and L. fasciculatum (1). 
Finally, eight plants exhibited external features typical of the genus 
Lithothamnion; however, none of these plants exhibited diagnostic 
characters necessary for their identification at the species level. 
Nevertheless, based on some morphological distinctions, they were 
partitioned into four morphospecies temporarily labeled as 
Lithothamnion spl (3 collections), Lithothamnion sp2 (1), 
Lithothamnion sp3 (1), and Lithothamnion sp4 (3). 

DNA extraction, PCR amplification and sequencing 

A subsample for DNA extraction was obtained by grinding a 
portion of the living surface of each specimen after avoiding areas 
with epiphytes, animal structures, and/ or damaged tissue. Special 
cautions were taken with the neotype of P. calcareum because this 
specimen has been archived in BM since 1983 [35]. To avoid 
contamination, this archival specimen was processed (DNA 
extraction and PCR amplification) individually with fresh batches 
of reactants on a separate date after carefully cleaning the 
laboratory. To increase the possibility of detecting contamination, 
several genes were amplified for this specimen on the same date 
(SSU, rhch, psbA, COI-5P) running negative controls in parallel 
for each gene; none of the chromatograms showed evidence of 
background signal and all negative controls were clean [37]. The 
holotype of M. sphaericum in SANT was collected shordy before 
the present study (October 2008) and processed alongside 
topotype material of the same species. Attempts to acquire 
sequence data from type material of other species included L. 
corallioides, L. fornicatum, N. hrassica-florida, and S. fruticulosa 
but proved unsuccessful. 

DNA was extracted with the DNeasy Blood & Tissue Kit Spin- 
Column Protocol (Qiagen) following manufacture's recommenda- 
tions. Two gene fragments were amplified: (i) a fragment of 664 bp 
of the standard DNA barcode (the 5' end of the mitochondrial 
gene cytochrome oxidase I, COI-5P) with primers GazFl and 
GazRl from [14], and GCorR3 (5'TGATTYTTYGGA- 
CATCCTGA3'), and (ii) a fragment of 892 bp of the plastidial 
gene photosystem II reaction center protein Dl (pshA) with 
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Figure 1. Primary Species Hypothesis (PSHs) delineated witli tiie COI-5P gene. Bayesian gene tree with posterior probabilities (>0.9) next 
to each node. Branch tips are the 29 haplotypes detected in the study. Vertical thick lines indicate PSHs delineated with ABGD and GMYC methods; 
numbers next to the vertical lines are PSH codes; letter A-D next to some nodes indicate major lineages. For ABGD, partitions for the more inclusive 
(M) and less inclusive (L) results are shown. GMYC partitions include the single-threshold (S) and multiple-threshold (M) variants of the method. Grey 
thick lines indicate discrepancies between partitions. 
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primers pshA-Fl and psh A-K2 from [38]. PGR reactants were 
prepared in a laminar flow hood and PCRs were performed in 
25 |xL containing 2 |J,L of DNA template, 2.5 |lL of 1 x PGR 
buffer, 2.5 mM MgG12, 0.192 mM dNTPs, 0.1 nM of each 
primer, and 1.2 U of Taq DNA Polymerase (Sigma- Aldrich) in a 
Biometra TProfesional Basic thermocycler following [39]. Ampli- 
fication success was evaluated by electrophoresis. After removing 
the excess of primers and nucleotides with shrimp alkaline 
phosphatase and exonuclease I enzymes, PGR products were 
bidirectionaUy sequenced at Macrogen facilities (http:/ /www. 
macrogen.com). All sequences are publically available in BOLD 
and GenBank (see Table SI for BOLD IDs and GenBank 
accession numbers). 

Data analyses 

Sequences were aligned and edited using the program Geneious 
5.6.6. As we aimed to delimit species based on sequence data 
rather than to assess their phylogenetic relationships, we chose not 
to run maximum-likelihood or maximum parsimony analyses. 
Instead, GOI-5P sequences were partitioned into a set of PSHs 
using two bioinformatics tools: ABGD [30] and GMYC [29,31]. 
For ABGC, genetic distances between specimens were calculated 
using the Kimura two parameters (K2P) model, a standard metric 
in DNA barcoding studies. ABGD was remotely run at http:// 
wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html using default 
values except for the relative gap width {X) which was set to 10 
to avoid the capture of small local gaps. For the GMYC method, 
duplicate haplotypes were removed from the alignment using 
DnaSP 5.10.01 [40]. Since the GMYC method is based on 
branching rates, branch lengths were estimated under a relaxed 



log-normal clock with the Bayesian analysis implemented in 
BEAST 1.7.4. Following [31], BEAST was run using a coalescent 
(constant population size) prior and the best-fitting model 
identified by jModelTest (HKY-nG with G = 0.153) [41,42]; the 
parameters for the substitution model (substitution rate, rate 
heterogeneity, and base frequencies) were unlinked across 
positions. MGMC chains were run for 20 million generations 
with a 10% burnin (determined by visual inspection of MCMG 
progression). After termination, the MCMG output was analyzed 
with TreeAnnotator 1.7.4 using all trees after the burnin, a 
posterior probability limit of 0.5, targeting the maximum clade 
credibility tree, and keeping the target node heights. Both the 
single-threshold and the multiple-threshold versions of the GMYC 
model [29,31] were optimized onto the output tree with the help 
of the SPLITS v. 1.0-19 package for R. AlC-based support values 
for the GMYC clusters were calculated following [32]. BEAST 
and TreeAnnotator were also employed to reconstruct a 
phylogeny for the psbA gene with the same options used for the 
COI-5P sequences but a different best-fitting model (GTRh-G with 
G = 0.175). 

Results 

Primary Species Hypothesis delineation based on C0I-5P 
sequence data 

The 224 collections of maerl were sequenced for a 664 bp 
fragment of the barcoding COI-5P gene; 29 unique haplotypes 
were found with 227 variable sites. Genetic pairwise K2P distances 
ranged from 0 to 0.21 while the shape of the pairwise distance 
distribution was clearly bimodal with two conspicuous peaks at 
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pairwise distances <0.01 and 0.14-0.16 separated by a rough gap 
of very low frequencies. The number of PSHs delineated with the 
ABGD method varied with the maximum prior distance (P) used 
in the analysis. Extreme prior thresholds led to uninformative 
partitions where either each haplotype was delimited as a different 
species or all haplotypes were included into a single PSH. 
Intermediate values of P led to partitions with 9 (P = 0.013), 13 
(P = 0.008) and 14 PSHs (P = 0.0017 to 0.005). Partitions with 9 
and 14 PSHs are detailed in Fig. 1. 

GMYC was applied to a phylogenetic tree reconstructed with a 
relaxed lognormal clock. Effective sample size for each statistic of 
the tree always was >500 and the MCMC converged to a 
stationary distribution. The likelihood of the null model 
{Lq — 170.19) was significantly lower than the maximum likelihood 
of the single-threshold version of GMYC model (isingie= 177.97, 
P-value = 0.0004). According to the latter, the transition from 
speciation to coalescent occurred at a depth of 0.0014 substitutions 
per site and resulted in a partition with 13 PSHs (confidence 
intervEil 4—14): 7 distinct clusters plus 6 singletons (Fig. 1). The 
likelihood of the multiple-threshold version of the model 
(innjitipie^ 178.36) also was significandy higher than that of the 
nuU model (P-value = 0.00028). This version detected a second 
threshold for the speciation-coalescent transition towards the tips 
of the tree at an extremely shallow depth of only 0.00026 
substitutions per site. With tiiis new threshold, the analysis 
delimited 15 PSHs (confidence interval 4—15): 7 clusters plus 8 
singletons. The mean support value across GMYC clusters was 
similar in the single-threshold (0.73 ±0.1 38) and in the multiple- 
threshold (0.74±0.237) methods. Nonetheless, the two new 
clusters delimited by the multiple-threshold algorithm had very 
litde support (<0.45). 

Although based on entirely different criteria, the partitions 
delineated by ABGD and GYMC were notably congruent. The 
less inclusive partition obtained by ABGD (14 PSHs) was nearly 
identical to the one produced by the single-threshold version of the 
GMYC model (13 PSHs). The only discrepancy involved 
haplotype Hap_24 (a specimen from Svalbard Archipelago) which 
was resolved as a singleton by ABGD while the GMYC model 
clustered it with other collections from Svalbard and Scandinavia 
(Hap_7). In comparison, the more inclusive partition of ABGD (9 
PSHs) and the multiple-threshold GMYC (15 PSHs) showed more 
discrepancies. However, their conflicting PSHs (greyed in Fig. 1) 
seemed biologically implausible because (i) the more inclusive 
hypothesis of ABGD clustered groups of haplotypes separated by 
average K2P distances as large as 0.064 (PSHs 4 to 7 in Fig. 1) or 
0.081 (PSHs 10 to 12), and (ii) the conflicting PSHs in the multiple- 
threshold version of GMYC were separated by distances as small 
as 0.002-0.005 (sequence divergence between and within PSHs is 
shown in Table S2a). 

C0I-5P and psbA phylogenies and secondary species 
hypothesis 

The phylogeny inferred from the mitochondrial COT5P gene 
was well resolved (15 nodes with posterior >0.95 out of 26 nodes) 
and four major lineages could be distinguished (Fig. 1). Regardless 
of the approach used for species delineation, all the PSHs that 
included more than one COI-5P haplotype always coincided with 
clades with high statistical support (posterior >0.9) with the 
exception of PSH 7. 

The chloroplastic psbA gene was sequenced for fifteen of the 
sixteen PSHs defined with COI-5P data. A fragment of 892 bp 
generated for 29 specimens produced 15 haplotypes. The 
phylogeny inferred from psbA. data was remarkably congruent 
with the one inferred from the COI-5P gene (Fig. 2). Again, four 



major lineages could b(- n^cognized, one of them separated from 
the others at an earUer time. Eleven PSHs were characterized by 
unique pshh haplotypes. The remaining four PSHs shared pshA 
haplotypes by pairs (PSH 4-5, PSH 6-7); these pairs corresponded 
to those cases where the delineation produced by ABGD was in 
conflict with the solution of the GMYC method. Most psbA 
haplotypes exhibited pairwise distances within a range of 10 
(equivalent to 98.8% similarity) to 113 differences (86.6%). StiU, a 
few PSHs were characterized by psbA sequences separated by 
distances as smaU as 2 (PSH 13 vs. PSH 14, 99.8% similarity), 3 
(PSH 15 vs. PSH 16, 99.6%) or 5 point mutations (PSH 10 vs. 
PSH 11, 99.4%). Only two PSHs (PSH 3 and PSH 4) produced 
more than one (two) psbA haplotypes that were separated by a 
single mutation (99.9% similarity). 

Given the consistency between the two phylogenies, any PSH 
drawn using ABGD and GA4YC that was either monophyletic for 
pshA or had unique pshA haplotypes was proposed as a SSH. As a 
result, 16 initial PSHs were converted to 13 SSHs (Fig. 2; see also 
Table S2b for average COT5P sequence divergence between and 
within SSHs). Two pairs of PSHs alternatively recognized as either 
a single PSH or as two different PSHs by ABGD and GMYC were 
turned in a single SSH each (SSH 4+5, SSH 6+7). Unfortunately, 
we did not manage to obtain a pshA secjuc-ncc- for PSH 9, a PSH 
supported by only one of the partitions derived from the GMYC 
model. Since its COI-5P sequence was very close to the haplotypes 
found in PSH 8, we opted for a conservative inclusive approach 
and considered these two PSHs as a single SSH (8+9). The final 
partition into SSHs matched the delineation obtained with the 
single-threshold alternative of GMYC applied to COT5P data 
only and was nearly identical to the less inclusive partition 
generated by ABGD. 

Matches in public data bases 

Based on the literature and on the magnitude of intra-SSH 
variability found in our study, we used ad-hoc cutoff values (>98% 
identity for COI-5P, >99% for psbA) to determine which 
GenBank searches had returned hits for potential conspecifics. 
Only 9 out of our 29 COI-5P haplotypes (4 SSHs) and 3 out of the 
15 psbA haplotypes (2 SSHs) resulted in a relevant match in either 
GenBank or BOLD (Table S3). Altogether, we obtained hits for 4 
out of our 13 SSHs: SSH4+5, SSH6+7, SSH 12, and SSH 16. Only 
SSH4+5 resulted in a match to an identified species, L. glaciale: 
our COT5P sequences were 98.6-99.9% similar to, and shared 
the same Barcode Index Number (BOLD:AAA6958), 39 acces- 
sions uploaded to BOLD for plants collected in Northeast USA 
and Canada. Also from OSPAR region I, SSH6+7 was conspecific 
with plants (BIN BOLD: ABA9580) from the Pacific (British 
Columbia) which, according to pictures logged in BOLD, have a 
branched morphology typical of maerl-forming plants. Finally, 
SSH 12 and SSH 16 had conspecific matches in GenBank with 
specimens from Brittany, France, which were only identified as 
CoraUinales. 

Concordance with morphological identification 

The total number of species identified based on their 
morphological features (10) was close to the number of SSHs 
(13) delimited with molecular data. However, among the 11 SSH 
with more than one specimen, only SSHs 1, 3, 11, and 13 were 
consistently assigned to a single morphospecies (Fig. 2). Many 
morphospecies contained collections from two, five or even six 
distinct SSHs. The only exceptions were M. sphaericum, a maerl 
species with a distinctive spherical morpholog)', and the three 
collections assigned to morphospecies Lilholhamnion sp2 that 
clustered under SSH 2. 
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PSH SSH HYPOTHESIS BASED ON MORPHOLOGICAL FEATURES DISTRIBUTION 

IV 



I 
I 



I I Ulholhamnion airalliodes 
I I Lilhothmmon gkciale 
I I Lilholhamnm sp\ 
I I Lilhothmmon sp2 
Lilhollinmnion spj 
[ I Lilholhomiiion st>4 



^1 Mesophyllm sphamivm 
Lilhophllm fasciaihiiiin 

□ UthofMlm 

Phynialolilhon calcoKum 



TAXONOMIC CONCLUSIONS 
MesopMhm ^phmiam ** 

Mesophyllm sp2 



Lahollwjiiiiioit glcidale * 

Litholhamm sp2 
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Phyimfoljrlwii 
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Phymilolilhoii 5p3 

Phymlolllhon sp4 
Litliophyilwii fasncniutiim 
Llrhophylli/iu Jeii/almi 



Figure 2. Secondary Species Hypothesis (SSHs) corroborated with the p5i&A gene. Bayesian gene tree with posterior probabilities (>0.9) 
next to each node. Branch tips are the 16 psbA haplotypes detected in the study. Numbers at the tip of the branches are Primary Species Hypothesis 
(PSH) codes (see Fig. 1) while vertical thick lines delineate SSHs; letters A-D next to some nodes indicate major lineages. Stacked horizontal bars next 
to the gene tree indicate the morphospecies identified at the onset of the study; numbers within bars are the actual number of specimens recorded 
for each morphospecies. The distribution of each SSH across OSPAR regions and its taxonomic identity is also provided. * topotype specimens 
sequenced, ** type specimen sequenced. 
doi:1 0.1 371 /journal.pone.01 04073.g002 



Attribution of available species names 

Six SSHs could be assigned to a species name using a body of 
proofs. In two cases, name assignment rested on comparisons with 
molecular data obtained from type material. SSH 12 was identical 
to COI-5P sequences obtained from neotype material of P. 
calcareum from BM. Likewise, our collections of SSH 1 included 
the holotype of M. sphaericum. 

For the second most widespread and common species in our 
study (SSH 8+9), we tentatively attributed the species name L. 
corallioides because the latter, together with P. calcareum, is 
typically regarded as a common component of maerl beds in 
Atlantic Europe. Furthermore, samples from the type locality of L. 
corallioides (Rade de Brest, Finistere, France) [34] were resolved in 
the SSH 8+9. 

In light of the morphological traits observed by SEM together 
with the existence of previous records from the same area, we 
temporarily attributed the names L. fasciculatum and L. dentatum 
to SSH 15 and SSH 16. In doing so, we used two names currently 
available in the literature for the European Atlantic while 
acknowledging that they belong to entities that need revision. A 
reassessment of the lectotype of L. fasciculatum has revealed that 



the epithet fasciculatum was misapplied to Atlantic plants 
belonging to the genus Lithophyllum [43]. Likewise, it seems 
unlikely that the coraUine algae identified as L. dentatum in the 
Atlantic and their Mediterranean counterparts may be conspecif- 
ics [44] . Indeed, the Atlantic plants of L. dentatum were previously 
considered a form of L. incrustans Philippi [45] . 

Finally, we attributed the species name L. glaciale to SSH 4+5 
based on the result of the BOLD identification engine. The name 
used in BOLD has not been confirmed by matching to sequences 
of type material (see [46]) and should be used with caution. 
Nonetheless, we also recorded SSH 4+5 in Spitsbergen Island 
(Svalbard Archipelago), the type locality of L. glaciale and where 
this coralline is reported to be common along the west and north 
coasts of the island (see [36] and references therein). We did not 
dare to link other SSHs to available species names and, therefore, 
seven SSHs were left without a binomial name. Nevertheless, their 
generic affiliations were evident based on their morphological 
traits and the phylogenetic relationships inferred from our COI-5P 
and psbA data, and we temporarily named them as Lithothamnion 
sp.2 (SSH 6+7), Phymatolithon sp.l (SSH 10), Phymatolithon sp.2 
(SSH 11), Phymatolithon sp.3 (SSH 13), Phymatolithon sp.4 (SSH 
14), Mesophyllum sp.l (SSH 2), and Mesophyllum sp.2 (SSH 3). 
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Figure 3. Morphological variability and distribution of the Secondary Species Hypotheses (SSHs) along the OSPAR area. The table 
indicates the number of specimens assigned to each SSH per OSPAR region and country. Vertical lines delimitate four major lineages revealed by 
both C0I-5P and psb/K phylogenies; tentative names are provided for the clades (genus level) based on sequence information from type/topotype 
material and on the occurrence of conspicuous morphological features. Lithoth. = Lithothamnion, Phymat. = Phymatolithon, Lithoph. = 
Lithophyllum, IVIesoph. = Mesophyllum. Scale divisions in the photographs are mm. 
doi:1 0.1 371/journal.pone.01 04073.g003 



Geographical distribution 

A majority of SSHs were restricted to one (6 SSHs) or two (4 
SSHs) sampling areas (Fig. 3). The remaining three SSHs showed 
wider, largely continuous distributions. Lithothamnion sp.2 (SSH 
6+7; for equivalence between SSHs and species names see Fig. 2) 
seemed confined to high latitudes in OSPAR regions I and II 
where it showed minimal overlap with other maerl-forming plants 
except the phylogenetically close L. glaciale (SSH 4+5). South- 
ward, Lithothamnion sp.2 was replaced by two species, L. 
corallioides (SSH 8+9) and P. calcareum (SSH 12), with wide 
rang(;s that reached northwest Spain. None of the species found in 
Macaronesia (Madeira and Canary Islands) was detected in the 
OSPAR area and the other way around. 

Each major lineage resolved in our phylogeny had different 
ranges and limits. Lithothamnion (lineage A, Fig. 3) was detected 
at mid to high latitudes in the OSPAR area, with a southern edge 
in northwest Spain. Phymatolithon (Lineage B) went unrecorded in 
OSPAR region I; instead, it reached the warmer coasts of southern 
Portugal and Macaronesia, albeit with a replacement of species. 
Finally, Lithophyllum (lineage C) was restricted to the British Isles 
and north of France whereas Mesophyllum (lineage D) was 
confined to southern latitudes (Spain, Portugal and Macaronesia). 

Discussion 

Delimitation and identification of species of coralline 
algae forming maerl 

Our spatially comprehensive sampling likely provides a 
thorough picture of the alpha diversity of maerl along European 
Adantic coasts. Our analyses of both COI-5P and psbA sequence 
data have delineated 13 SSHs, a number comparable to the 16 
species reported for the OSPAR regions and Macaronesia (see 
Table 1 for references). Nevertheless, linking available binomial 
names to the SSHs uncovered in our study was a delicate task. In 
addition to various analyses of the sequence data, additional 
evidence (morpho-anatomical observations, previously known 
species distribution, molecular data from type and topotype 
specimens) was required to guide our decisions at the time of 
attributing names. Following this approach, we managed to name 
almost half of the species detected in our study with acceptable 
confidence. However, it is likely that most, if not all, of the species 
that we left unnamed in this study may have already been 
described elsewhere. 

We did not dare to identify SSH 6+7 (Lithothamnion sp. 2) to 
species level. Its confinement to OSPAR region I and to 
northernmost sites of region II suggests that L. tophiforme could 
be a plausible name but our sequences did not match two 
collections from New Foundland, Canada lodged in BOLD as L. 
tophiforme. However, the identity of these collections in BOLD 
has not been confirmed with sequence data from type material and 
should be considered with caution. Lithothamion tophiforme is 
mainly reported as an Arctic species that, in European waters, is 
confined to very high latitudes [47]. Somewhat unexpectedly, 
however, L. tophiforme went unrecorded in a detailed recent 
investigation of the northernmost maerl communities currently 
known, discovered in 2006 at 80°3rN in the Svalbard Archipel- 
ago [36]. Instead, these communities were dominated by L. 



glaciale, the only species that we found in our collections from 
Svalbard. Interestingly, our BOLD searches revealed that 
Lithothamnion sp. 2 also occurs in the North Pacific (British 
Columbia). A comparable circumpolar distribution has been 
reported for L. glaciale (see references in [48]). 

Our results show that Phymatolithon sp.3 (SSH 13) is a major, 
even dominant, component of maerl beds in Spain and Portugal. 
Indeed, a recent quantitative study with DNA barcodes demon- 
strates that the widespread belief that L. corallioides and P. 
calcareum are the major builders of maerl in the temperate 
European Atlantic does not hold for the Iberian Peninsula. 
Instead, they are gradually replaced by Phymatolithon sp.3 in 
Galicia (NW Spain) to become extremely rare in S Portugal [49]. 
Despite our efforts, we have been unable to resolve the identity of 
this species beyond generic level. The examination by SEM 
revealed traits also found in Phymatolithon lamii (e.g. sunken, 
rimless conceptacles), a common coralline throughout the British 
Isles, northern Spain, France, Norway, Iceland and eastern North 
America [50]. It has also been reported from the western North 
Pacific Ocean [51,52] and, more recently, from the Mediterra- 
nean where it might be an alien species [53]. In addition, one of 
the co-authors (V.P., unpubfished results) recently sequenced a 
600 bp long fragment of the pshh gene from the type specimen of 
P. lamii (in PC) that reveals a low-moderate divergence with our 
Phymatolithon sp.3. However, P. lamii has always been described 
as encrusting thaUi, and there is no previous record of its 
occurrence as maerl. Hence, further sequence data from type 
material will be required to assign a species name to SSH 13. 
Lastly, our results clearly indicate that the maerl-forming algae 
that colonize Macaronesia deserve further study with appropriate 
sampling design and molecular tools. 

To our knowledge, this is the first barcoding study focused on 
maerl-forming algae. Previous DNA barcode studies on coralline 
algae mostly focused on geniculate forms [13,22,46,54] or were 
intended to resolve infra-ordinal phylogenetic relationships among 
the CoraUinales [55]. Nevertheless, Bittner et at. [56] sequenced 
geniculated and non-geniculated coralline algae, mainly from 
South Pacific Islands, for pshA and COI-5P. and used ABGD and 
GYMC to delineate spc'cies. As the authors found very divergent 
numbers of 'genetic spetdes' depending on the criteria and the 
marker, they concluded that DNA-barcoding was non-accurate for 
assessing the species diversity in this group (but see [33] for 
alternative explanations when GMYC has a highly divergent 
outcome). Contrarily, we are in favor of an integrative systematic 
approach to investigate the diversity of maerl-forming red algae. 
We propose the use of the mitochondrial COI-5P barcode as the 
first marker and the plastidial psbA gene as a secondary marker 
along with other lines of evidence. In this regard, we follow other 
authors that already noted the intrinsic limitations of delimiting 
species from single-locus studies and advocate the incorporation of 
multiple lines of evidence (biogeographical, biological, additional 
gene sequences) in this studies [10,32,33]. 

Despite considerable efforts to identify specimens collected in 
this study based on their morpho-anatomical characters, sequence 
data were incomparably more efficient. In fact, our Fig. 2 shows 
that even maerl assigned to different morphogenera turned out to 
belong to the same molecular entity. Initially, this considerable 
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discrepancy between the morphological identification and the 
molecular-based delimitation might seem astonishing. However, it 
is just another example of the considerable challenge of 
discriminating species in a group that largely lacks diagnostic 
features but shows considerable morphological plasticity and 
convergence. Species delimitation with molecular markers is 
known to perform well with deeply divergent taxa such as those 
found in our study. However, we stiU uncovered a few closely 
related species that deserve further studies. It is increasingly 
acknowledged that a multi-locus approach can define reliable 
species boundaries in those cases [10,29,31,57,58]. Even for the 
more closely related SSHs, the genealogical concordance observed 
between loci located in different genomes lends additional support 
to our partition. Finally, the congruence between partitions 
recovered with analyses based on radically different criteria 
(coalescence vs. distribution of pairwise differences between 
sequences) indicates that there is a strong signal in the COI-5P 
data set. If any, the only shortcoming encountered in the course of 
that study was the erratic amplification of the COI-5P fragment 
that failed on 1/3 of the specimens at the first attempt, although 
repeated amplifications and/or re-extracting DNA from the same 
individual often solved this issue. The remarkable coincidence 
between the patterns revealed by our pshK and COI-5P sequence 
data indicates that the plastidial marker, which is easier to amplify 
than COI-5P, could be a useful alternative for the identification of 
maerl-forming species (but see [54] for a intrageneric study where 
the gene trees produced by COI-5P and by another plastid- 
encoded gene, rficL, were not always congruent). Compared to the 
standard DNA barcode, jf),s/;A has lower resolution. However, most 
of the maerl-forming taxa studied in the OSPAR region are deeply 
divergent. The only exceptions were two pairs of closely related 
species with low levels of interspecific divergence (2-3 nucleotide 
substitutions) that could be mistaken for intraspecific variability (L. 
fasciculatum vs. L. dentatum; Phymatolithon sp.3 vs. Phymatoli- 
thon sp.4). 

Species distribution and implication for future prospects 

Our study reveals that two species of coralline algae are the 
main constituent of maerl beds in temperate European Atlantic: L. 
corallioides and P. calcareum. Another Phymatolithon (SSH13) 
replaces them in the south while the cold OSPAR region I seems 
dominated by two species of Lithothamnion [L. glaciale and 
Lithothamnion sp.2-SSH 6+7). The remaining species unraveled in 
our study are either infrequent and/or confined in space. The 
gradual replacement of species with latitude (Fig. 3) is consistent 
with the patterns of maerl distribution reported in the literature 
(see Table 1 and references therein). The distribution of coralline 
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